
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.1
 *  @date    Tue Sep 17 17:12:07 EDT 2013
 *  @see     LICENSE (MIT style license file).
 */

package scalation.linalgebra

import java.lang.Math.copySign

import math.{abs, sqrt}

import scalation.linalgebra.MatrixD.eye

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `Givens` objects has methods for determinng values 'c = cos(theta)' and
 *  's = sin(theta) for Givens rotation matrices as well as methods for applying
 *  Givens rotations.
 */
object Givens
{
    type CosSin = Tuple2 [Double, Double]    // type for (cosine, sine)

    private val DEBUG = false                // debug flag

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create the values for a Givens 2-by-2 rotation matrix.  Given scalars
     *  'y' and 'z', efficiently compute 'c = cos(theta)' and 's = sin(theta)'
     *  that can be used to form the rotation matrix.
     *  @see Algorithm 5.1.3 in Matrix Computation.
     *  @param y  the first scalar
     *  @param z  the second scalar
     */
    def givens (y: Double, z: Double): CosSin =
    {
        if (z == 0.0) {
           (1.0, 0.0)
        } else if (abs (z) > abs (y)) {
           val t = y / z
           val s = -copySign (1.0 / sqrt (1.0 + t * t), z)
           (-s * t, s)
        } else {
           val t = z / y
           val c = copySign (1.0 / sqrt (1.0 + t * t), y)
           (c, -c * t)
        } // if
    } // givens

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Efficiently perform a Givens row update: 'a = g (i, k, theta).t * a'.
     *  The update just affects two row.  FIX
     *  @param a   the matrix to update
     *  @param i   the first row ??
     *  @param k   the second row ??
     *  @param cs  the (cosine, sine) of theta
     */
    def givensRowUpdate (a: MatrixD, i: Int, k: Int, cs: CosSin)
    {
        for (j <- 0 until a.dim2) {
            val r1 = a(i, j)
            val r2 = a(k, j)
            a(i, j) = cs._1 * r1 - cs._2 * r2
            a(k, j) = cs._2 * r1 + cs._1 * r2
        } // for
    } // givensRowUpdate

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Efficiently perform a Givens column update: a = a * g (i, k, theta).
     *  The update just affects two columns.  FIX
     *  @param a   the matrix to update
     *  @param j   the first column ??
     *  @param k   the second column ??
     *  @param cs  the (cosine, sine) of theta
     */
    def givensColUpdate (a: MatrixD, j: Int, k: Int, cs: CosSin)
    {
        for (i <- 0 until a.dim1) {
            val r1 = a(i, j)
            val r2 = a(i, k)
            a(i, j) = cs._1 * r1 - cs._2 * r2
            a(k, j) = cs._2 * r1 + cs._1 * r2
        } // for
    } // givensColUpdate

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a Givens rotation matrix with angle 'theta = atan(s/c)'.
     *  A matrix is post-multiplied by the Given matrix to clear element (i, k).
     *  The 2-by-2 rotation is embedded in an identity matrix of dimension 'n'.
     *  @param i   the first diagonal position  (i, i)
     *  @param k   the second diagonal position (k, k)
     *  @param n   the dimension of the resulting rotation matrix
     *  @param cs  the (cosine, sine) of theta
     */
    def givensRo (i: Int, k: Int, n: Int, cs: CosSin): MatrixD =
    {
        val b = eye (n)
        b(i, i) =  cs._1; b(i, k) =  cs._2         // |  c   s  |
        b(k, i) = -cs._2; b(k, k) =  cs._1         // | -s   c  |
        if (DEBUG) println ("givensRo: b = " + b)
        b
    }  // givensRo

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a transposed Givens rotation matrix with angle 'theta = atan(s/c)'.
     *  A matrix is pre-multiplied by the Given matrix to clear element (k, i).
     *  The 2-by-2 rotation is embedded in an identity matrix of dimension 'n'.
     *  @param i   the first diagonal position  (i, i)
     *  @param k   the second diagonal position (k, k)
     *  @param n   the dimension of the resulting rotation matrix
     *  @param cs  the (cosine, sine) of theta
     */
    def givensRoT (i: Int, k: Int, n: Int, cs: CosSin): MatrixD =
    {
        val b = eye (n)
        b(i, i) = cs._1; b(i, k) = -cs._2           // |  c  -s  |
        b(k, i) = cs._2; b(k, k) =  cs._1           // |  s   c  |
        if (DEBUG) println ("givensRoT: b = " + b)
        b
    }  // givensRoT

} // Givens object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `GivensTest` object tests the `Givens` object by using two rotations to
 *  turn matrix 'a' into an upper triangular matrix, clearing positions (1, 0)
 *  and (2, 1).  A third rotation clears position (0, 1) but also puts a non-zero
 *  value at (1, 0).
 *  @see http://en.wikipedia.org/wiki/Givens_rotation
 */
object GivensTest extends App
{
    import Givens._

    var a = new MatrixD ((3, 3), 6.0, 5.0, 0.0,
                                 5.0, 1.0, 4.0,
                                 0.0, 4.0, 3.0)
    println ("a = " + a)

    val cs1 = givens (a(0, 0), a(1, 0))    // rotation 1 takes (1, 0) = 5 -> 0 
    val g1  = givensRoT (0, 1, 3, cs1)     // use givensRoT below the diagonal
    a = g1 * a                             // pre-multiply by Givens matrix g1
    println ("(c1, s1) = " + cs1)
    println ("rotation 1: a = " + a)

    val cs2 = givens (a(1, 1), a(2, 1))    // rotation 2 takes (2, 1) = 4 -> 0 
    val g2  = givensRoT (1, 2, 3, cs2)
    a = g2 * a                             // pre-multiply by Givens matrix g2
    println ("(c2, s2) = " + cs2)
    println ("rotation 2: a = " + a)

    val cs3 = givens (a(0, 0), a(0, 1))    // rotation 3 takes (0, 1) = 4.48 -> 0 
    val g3  = givensRo (0, 1, 3, cs3)      // use givensRo above the diagonal
    a = a * g3                             // post-multiply by Givens matrix g3
    println ("(c3, s3) = " + cs3)
    println ("rotation 3: a = " + a)

} // GivensTest object

